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ABSTRACT 

We carry out N-body simulations of several non-Gaussian structure formation models, 
including Peebles' isocurvature cold dark matter model, cosmic string models, and a 
model with primordial voids. We compare the evolution of the cluster mass function 
in these simulations with that predicted by a modified version of the Press-Schechter 
formalism. We find that the Press-Schechter formula can accurately fit the cluster 
evolution over a wide range of redshifts for all of the models considered, with typical 
errors in the mass function of less than 25%, considerably smaller than the amount 
by which predictions for different models may differ. This work demonstrates that the 
Press-Schechter formalism can be used to place strong model independent constraints 
on non-Gaussianity in the universe. 

Key words: Cosmic strings - cosmology: theory - large-scale structure of universe 
- galaxies: clusters: general 



1 INTRODUCTION 

Many cosmological models, including 'inflation', predict that 
the large scale structure of the universe arose from the grav- 
itational collapse of Gaussian primordial fluctuations. These 
models have a strong physical motivation (in inflation, the 
primordial perturbations are just quantum fluctuations), re- 
quire few free parameters (the two-point correlation func- 
tion provides a complete statistical description of the field), 
and have many observable properties which can be modeled 
with simple physics. Given these attractive features, Gaus- 
sian models are a natural starting point for the investigation 
of structure formation. 

A number of theories however, including 'defects' (Kib- 
ble 1976, Vilenkin & Shellard 1994) and certain exotic forms 
of inflation (Peebles 1983, 1997, 1998a,b; La 1991; Amen- 
dola & Occhionero 1991; Amendola & Borgani 1993), predict 
that the primordial fluctuations were not Gaussian. These 
theories possess many of the attractive features of Gaus- 
sian models: defect models, for instance, are well motivated 
and require few free parameters. Unlike the Gaussian case, 
however, making robust predictions in non-Gaussian mod- 
els is often very difficult. In defect theories, this is because 
the physical processes which must be modeled are complex 
and highly non-linear. More generally, the parameter space 
available to non-Gaussian models is infinitely larger, because 
specifying the statistics of the primordial fluctuations re- 
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quires not just the two-point correlation function, but all 
higher n-point correlation functions as well. 

Previous studies of non-Gaussian structure formation 
(e.g., Weinberg & Cole 1992; Borgani et al. 1994; Park, 
Spergel & Turok 1991) have tended to concentrate on spe- 
ciflc models, making predictions for observable properties 
in each case. Due to the infinite range of possible parame- 
ters, and the difficulties of making robust predictions even 
within well specified models, we turn here to a different ap- 
proach: We start from a property of non- Gaussian models 
which can be fully quantified by a small number of parame- 
ters, the probability distribution function (PDF) of the den- 
sity field over the range of scales relevant to galaxy clus- 
ter formation. While knowledge of the PDF does not fully 
specify the non-Gaussian nature of the distribution (which 
would require knowledge of all the higher n-point correla- 
tion functions), we show that it provides us with sufficient 
information to make robust predictions for a very interest- 
ing observable quantity: the redshift evolution of the cluster 
number abundance. We do this by demonstrating that the 
cluster evolution in non-Gaussian N-body simulations can 
be accurately described by a version of the Press-Schechter 
formalism (Press & Schechter 1974) modified to allow for 
non-Gaussianity (Chiu, Ostriker & Strauss 1997), which re- 
quires only the PDF as input. The result of this work is a 
powerful tool with which we can use observations of galaxy 
clusters to constrain the primordial PDF in the universe, 
without reference to any specific non-Gaussian models or 
the uncertainties which may be inherent to them. The work 
presented here validates a number of previous studies which 
have used the non-Gaussian Press-Schechter formalism, but 
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have not tested it (Chiu et al. 1998, Robinson, Gawiser & 
Silk 1998, van de Bruck 1998). Other work (Robinson, Ga- 
wiser & Silk in preparation) discusses in detail the strong 
model-independent constraints we can place using this tool. 

In §|2| we describe the PS formalism, and we extend 
the derivation in the non-Gaussian case to allow for a scale 
dependent probability distribution function. Section ^ de- 
scribes our non-Gaussian models and the realization of the 
initial density fields in each case. In §^ we give details of our 
N-body simulations. Finally, in we present our results and 
in §^ we discuss our conclusions. 



2 MODIFIED PRESS-SCHECHTER 
FORMALISM 

The Press- Schechter (PS) formalism (Press & Schechter 
1974) is a simple semi-analytic tool for predicting the ex- 
pected number density of clusters in an evolved density field, 
given the linear initial conditions. Extensive work has shown 
that this formalism can give a good fit to the cluster abun- 
dance observed in N-body simulations of Gaussian struc- 
ture formation models. Though the PS formula has mainly 
been applied to Gaussian fiuctuations, there is a straightfor- 
ward way to generalize it to the non-Gaussian case: Here, 
the statistics of the underlying density field are quantified 
by specifying its PDF Pr{5), where pii{5)dS is the prob- 
ability that the primordial field at a given point in space 
has an overdensity between S and 5 + dS after a top hat 
smoothing of scale R. We can rewrite pr{S) in terms of 
a rescaled PDF PR{y) which has zero mean and rms one: 
pii{S) = Pr(& I (Jr) I (JR, where or is the rms overdensity 
in spheres of radius R. For many non-Gaussian models it 
should be a reasonable assumption that only the rms value 
of (5 varies as a function of scale (for instance, in defect 
theories where the fluctuations are layed down in a scale 
invariant manner), so that Pfl(y) is in fact independent of 
R. The Press- Schechter formula for this case has been de- 
rived by Chiu et al. (1998) and used to make predictions for 
a number of non-Gaussian scenarios (Robinson et al. 1998; 
van de Bruck 1998). 

In this work, we rederive the formula for the more gen- 
eral case that the form of the PDF does depend on scale, 
since we will see later that the assumption of scale indepen- 
dence is valid for some but not all of the models we consider. 
The Press-Schechter formalism starts from the assumption 
that clusters with mass greater than some value M at red- 
shift z form in regions of the primordial density field whose 
overdensity smoothed on scale Rm and extrapolated using 
linear gravity to redshift z is greater than some critical value 
(5c. The scale Rm is chosen to be the radius of a sphere con- 
taining mass M in the primordial homogeneous universe, 
which satisfies 

M = YPmiiL, (1) 

where pm is the comoving matter density of the universe. 
The critical threshold 5c is taken from an analytic solution 
for the collapse of a spherically symmetric overdensity, which 
gives 5c = 1.69 for a critical universe (for a derivation of the 
weak cosmological dependence, see Lacey & Cole 1993, Eke 
et al. 1996). 

From this starting point, it is possible to derive the 



cluster mass function. The probability VyR of a point in 
space forming part of a cluster with Lagrangian radius larger 
than R is equal to the probability of the density field, after 
smoothing on scale R, having an overdensity larger than the 
critical overdensity 5c. That is. 



V>R = 



PR{5)d5. 



(2) 



To obtain the probability Vr dR of a point in space forming 
part of a cluster with Lagrangian radius between R and 
R + dR, we differentiate the above expression with respect 
to R and take the absolute value, 

- (3) 



PRdR- 



dR 



V: 



>R 



dR. 



We can obtain the number density of such clusters n{R) dR 
by dividing by the cluster Lagrangian volume. 



n{R) dR - 



3/ 
AtvRS 



dR 



V. 



>R 



dR, 



(4) 



where we have also multiplied by a correction factor /, whose 
value we will fix by ensuring that final mass function ac- 
counts for the entire mass of the universe. The integrated 
mass function 7V> m is then given by 



3/ 
47r 



dR 



(5) 



To fix the value of the correction factor /, we use the fol- 
lowing expression for the integrated mass density pm of the 
universe: 



Mr n{R) dR, 



(6) 



where Mr is the cluster mass corresponding to Lagrangian 
radius R. Substituting using equation |l| and ^, we find 



/=(P>fl|«=o) 



(7) 



which does indeed give / = 2 for Gaussian fluctuations, as 
originally proposed by Press & Schechter. 

We comment briefly on some computational issues: We 
can rewrite the mass function in a convenient form using 
integration by parts. 



iV>Af = — 



3/ 

47r 



i?3 



^V>R 



(8) 



For each model and for each redshift we compute and store 
the function VyR for a range of values of R. Each computa- 
tion of A''>jv_f (eq. ^ then requires only a single integration 
over these stored values. Finally, we note that the redshift 
dependence of the mass function arises through the depen- 
dence of VyR on redshift, which for conciseness we have not 
labelled explicitly here. 



3 NON-GAUSSIAN MODELS 

For the purposes of this work, we consider three classes of 
non-Gaussian model; Peebles' isocurvature cold dark mat- 
ter (ICDM) model, cosmic string models, and primordial 
void models. These models are all physically motivated, and 
each one represents an interesting structure formation sce- 
nario in its own right. However, each model has a number 
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of uncertain features. In the cases of ICDM or Voids, there 
is a wide range of possible model parameters to consider, 
and in the case of cosmic strings, there are extreme compu- 
tational difficulties involved in realistically evolving a string 
network. Here we will not be concerned with these uncertain- 
ties, since our goal is simply to pick a diverse set of concrete 
non-Gaussian models, and use these examples to test the 
applicability of the non-Gaussian PS formalism. To the ex- 
tent to which the formalism works, we will have a tool with 
which we can make specific predictions for cluster evolu- 
tion in a wide class of non-Gaussian models, taking as input 
only the initial PDF pr{5) (which is what enters into the PS 
calculation; see §^). We will then be able to use cluster evo- 
lution to constrain the initial PDF of the universe without 
referring to any specific non-Gaussian models, a possibility 
which will be extremely useful given the difficulty of making 
robust predictions in so many non-Gaussian scenarios. 

We now explain the motivation for each of the mod- 
els, and present the computational details of realizing initial 
density fields in each case. 

3.1 Peebles ICDM model 

The first non-Gaussian structure formation scenario we 
consider is the Peebles ICDM model (Peebles 1983, 1997, 
1998a, b). This is an inflationary model which gives rise to 
non-Gaussian isocurvature fluctuations in the matter distri- 
bution. The cold dark matter (CDM) in this model is a field 
<j}{x) with mass density 

p(x) (X 0(x)2 (9) 

which survives as a relic at the end of inflation. Pertur- 
bations in this CDM field grow from quantum fluctuations 
which are frozen during the accelerated expansion of the uni- 
verse. The perturbations in (^(x) are Gaussian, with power 
spectrum 

P0 cx fc'"* (10) 

over a broad range of scales, where the spectral index 
depends on the details of the model. The corresponding fluc- 
tuations in p are non-Gaussian, and it is easy to show that 
they have an initial power spectrum 

PpXk""", (11) 
with 

nip = 3 -f 2m^. (12) 

Due to the isocurvature nature of the fluctuations, the 
CDM density does not evolve until the matter era, and then 
only on scales which are inside the horizon. Those modes 
which have crossed the horizon before matter-radiation 
equality will all start to grow at that time, while those modes 
which enter the horizon later will be increasingly suppressed. 
For our investigation of cluster formation we are interested 
mainly in the former, smaller scales, so it will be sufficient 
to assume scale-independent linear growth of the initial fluc- 
tuations. 

We consider two values of the initial spectral index of 
the Gaussian field m^. First, we consider a model (which 
we denote ICDM-2.4) with m0 = —2.4, in a background 
cosmology with fim = 0.2 and Qa ~ 0.8. This set of param- 
eters has been suggested as an interesting starting point by 



Peebles (1998b) since it gives power spectra in reasonable 
agreement with large scale structure and CMB data. Sec- 
ond, we consider a model (ICDM-2.0) with = —2.0, in a 
background cosmology with flm = 1 and S7a — 0. We do not 
expect this model to give a particularly good fit to the large 
scale structure data, but we introduce it since the degree 
of non-Gaussianity is intermediate between that of (ICDM- 
2.4) and a Gaussian distribution, and it therefore allows us 
to test the Press-Schechter formalism over a wider range of 
models. 

For each ICDM model, we generate a 128'' realization 
of a Gaussian random field ^(x) with the appropriate power 
spectrum P{k) cx m^. We then compute p{x) = 4>{x)^, and 
the corresponding density perturbation 

S{^) = '-^, (13) 
where p is the mean density, given by 

JV 

and A*' is the number of points on the lattice. Finally, we 
rescale the density contrast so that the density field will 
have the correct value of erg when linearly evolved to z = 0. 

3.2 Cosmic Strings 

The second class of non- Gaussian model we consider is that 
of cosmic strings. Here, structure is seeded gravitationally by 
an evolving network of string-like 'defects', topological relics 
of a symmetry breaking phase transition in the very early 
universe (Kibble 1976; Vilenkin & Shellard 1994). Recent 
work suggests that cosmic strings in a critical universe do 
not give a viable theory of structure formation (Allen et al 
1997; Albrecht, Battye & Robinson 1997), although strings 
in a universe with a significant cosmological constant may be 
able to give a good fit to existing large scale structure data 
(Battye, Robinson & Albrecht 1998; Avehno et al 1998). For 
the purposes of this work, we shall not be concerned by the 
many uncertainties in the details of realistic string dynamics. 
Instead, we shall take a simple representative model and use 
it as a testing ground for our study of cluster evolution. 

Our string-seeded density models use the string simula- 
tions of Ferreira (1995). These simulations evolve a network 
of strings in flat space, where an exact solution of the equa- 
tions of motion exists (Smith & Vilenkin 1987). Although 
the string simulations are carried out in a flat space back- 
ground, the network evolution demonstrates many of the 
important properties which are seen in expanding universe 
simulations, the most important of which is the 'scaling' of 
the string density, ensuring that the induced density pertur- 
bations will be roughly scale invariant on the largest scales. 

Following Veeraraghavan & Stebbins (1990), the matter 
overdensity 5 induced by the network can be computed by 
convolving a component 9+ — Ooo + da of the string stress- 
energy with a suitable Green function, and integrating over 
time. That is, 

/""■ - 1 

5(k)=4^y dVr(k,r,')e+(k,V)Y-^^^-y^, (15) 

where rji and rjf are the initial and final times in the simula- 
tion, 5(k) is the Fourier transform of the induced overden- 
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sity, T2 is a transfer function taken from Veeraraghavan & 
Stebbins (1990), = /J dt/ a is the conformal time, a is the 
cosmological scale factor, and the last term is a 'compensa- 
tion factor', with compensation scale kc = 2-n /t), included to 
ensure artificially that the density evolution conserves stress 
energy on the largest scales. 

We carry out our flat space string simulations on a 128^ 
lattice, choosing the side of our box to be Mpc. We 

work out the conformal time 77 at each time-step by equat- 
ing it to the simulation horizon size (which is zero at the 
initial time and grows by one lattice spacing at each time- 
step). At each step, we compute the appropriate compo- 
nent S+(x) of the string stress energy tensor on the lattice, 
Fourier transform to obtain 6l+(k), and multiply by T(k, rf). 
The final density contrast is then given by the sum of indi- 
vidual contributions from each time-step, as in eq. ^ Since 
all significant perturbations are induced at a time when the 
fluctuations are still linear, it is sufficient to use linear grav- 
ity in this way to compute the initial density field, and only 
consider the effects of non-linear gravity in the subsequent 
evolution. We choose two representative string models. First, 
we consider a model with an fim = 1 hot dark matter (IfDM) 
background with one species of massive neutrino accounting 
for the entire mass of the universe, and a Hubble constant 
Ho = 100/i km s"^ Mpc"^ with h = 0.5. We denote this 
model 'Strings-HDM'. Second we consider a string network 
in an Q.m = 1 CDM background, also with h = 0.5, which 
we denote 'Strings-CDM'. 

3.3 Primordial Voids 

The final non-Gaussian model we consider is one where the 
matter distribution has an enhanced network of primordial 
voids. The physical motivation for such models comes from 
an infiationary theory with two scalar fields (La 1991): one 
field drives infiation, while at the same time the second field 
undergoes a first order phase transition. The bubbles nucle- 
ated in this phase transition give rise to under-densities in 
the CDM distribution, superimposed on Gaussian fiuctua- 
tions which are also generated during the inflationary epoch. 
The detailed properties of the bubble distribution, includ- 
ing their sizes, shapes and profiles, will depend on the exact 
nature of the infiationary model, but for some range of pa- 
rameters it is possible to produce bubbles which are large 
enough to have a significant impact on structure formation, 
yet still escape limits arising from current Cosmic Microwave 
Background (CMB) observations (Amendola & Occhionero 
1991; Amendola & Borgani 1993). 

A typical cosmologically interesting scenario would pro- 
duce a dense network of bubbles with radii of order 10-40/i~^ 
Mpc. For the purposes of this work, we generate a simple re- 
alization of such a model as follows: in a lOO/i"^ Mpc box, we 
generate a 128'^ realization of a Gaussian density field (5g(x), 
with a power spectrum Pp oc k~^'^. This spectral index is 
chosen so that the density field has a very similar power 
spectrum to that in the ICDM-2.4 model, thus ensuring that 
our simulations primarily compare the effect of varying the 
amount of non-Gaussianity in the fluctuations, rather than 
that of varying the spectral index. We then superimpose on 
this a distribution of Ny randomly located voids, where the 
voids are spherically symmetric top-hat under-densities with 
radius R = 12h~^ Mpc. The total density field is given by 



'5(x) = <5g(x)- 



JVv 

<5v^W^fl(x 

1=1 



Yi 



(16) 



where 5v is the amplitude of the void under-density, yi is 
the centre of the i"^ void and W_r(x) is the window function 
in real space of a spherical top-hat with radius R, that is: 



Wfl(x) = 







•ixj <i? 
. otherwise. 



(17) 



We dial the amplitude of the void under-density 5v to en- 
sure that the resulting matter field has a significant degree 
of non-Gaussianity, without completely swamping the Gaus- 
sian component of the distribution. We find Sy = O.lag.s is 
a suitable value, where CTg,8 is the rms value of the Gaus- 
sian component, after top-hat smoothing on a scale of 8h^^ 
Mpc. We choose Ny = 12, and generate the void centres yi 
by picking points at random, discarding any choices which 
would give rise to overlapping voids, and repeating until the 
number of acceptable centres equals Ny. Finally we rescale 
the field 5{x) to ensure that it has the chosen normalization 
as. 



4 DETAILS OF N-BODY RUNS 

In §^ we have explained the motivation for our non-Gaussian 
models, and set out in detail the processes involved in gen- 
erating density fields which realize these models. We now 
describe the techniques involved in transforming these linear 
density fields into initial conditions for N-body simulations. 
For each of the models in question we start from a density 
field 5(x) defined on a 128^ lattice. The N-body simulations 
require initial conditions in the form of positions and ve- 
locities of some number A^'p — rip particles, which we can 
generate from our density fields using the Zeldovich approx- 
imation (Zeldovich 1970). Specifically, we define a displace- 
ment field p(q) via 



V-p(q) = ~5o(q). 



(18) 



where the subscript zero denotes the density field extrapo- 
lated to the present day assuming linear growth. Starting 
from a regular grid of particles, the position x{z) of each 
particle at redshift z is given by 



x(z) =q + p(q)^. 



(19) 



where b{z) is the linear growth factor at redshift z, and q is 
the initial coordinate of the particle. The comoving velocity 
of that particle is given by 



v(q,2) = p(q) 



b{z) din 6 



(20) 



6(0) dt 

where t is the physical time. We can solve eq. |l^ for the 
displacement field efficiently using Fourier transforms, that 



p(k) = G(k)H/(k)5(k), 



(21) 



where p(k) is the Fourier transform of the displacement field 
p(x), G is the appropriate Green function, 



G(k) = -z-, 



(22) 
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and H^(k) = w{ki)w{k2)w{k3) is a smoothing function, 
whose form is chosen to optimize agreement between the 
power spectrum of the input density field (5(x) and that of 
the mass distribution represented by the displaced particles 
(which we can measure by resampling the particle distribu- 
tion as a new density field on a lattice). We find the choice 



w[k) — 



. kL/Up < 77 

. otherwise 



(23) 



gives good agreement. 

We evolve each density field using P^M simulations car- 
ried out on a GRAPE-3AF board attached to a Sun work- 
station. The GRAPE-3AF is a special purpose computer 
which quickly evaluates pairwise forces, allowing a cosmo- 
logical simulation with 64'' particles to be run in under a day. 
The code used for our simulations is based on the P3MG3A 
algorithm described by Brieu, Summers, & Ostriker (1995). 



Each simulation contains 



= nl 



64"^ particles with 

a particle mass of 1.05 X lO^^nmh'^ Mq. The effective force 
resolution is 0.128 times the mean inter-particle separa- 
tion, or 200/i"^ kpc for our lOO/i"^ Mpc box length. The 
initial redshifts zi for each model, chosen in each case by en- 
suring that < 0.15 at the Nyquist frequency, are given in 
Table |l|. For a given amplitude of fiuctuations, the criterion 
of no shell crossing generally forces us to adopt an earlier 
initial redshift in the non-Gaussian models than in the cor- 
responding Gaussian case. Outputs from the simulations are 
stored and analyzed for a number of redshifts between z = Zi 
and z — 0. 

For each output, we identify clusters in the simulations 
using the 'friends-of-friends' (EOF) algorithm. This algo- 
rithm defines clusters as those groups of particles which are 
connected by pairs whose separation is smaller than some 
linking parameter I, which is quoted in units of the mean 
inter-particle separation. We choose I — 0.2, so that the 
mean overdensity within the clusters identified is of order 
200. This choice of linking length has been shown to pick 
out clusters whose mass function is in good agreement with 
the PS prediction in the case of Gaussian fluctuations (Mo, 
Jing & White 1996). 



5 RESULTS 

We now compare the predictions of the modified PS formal- 
ism with the results of our simulations. In order to compute 
the PS prediction, we require as input for each model the 
amplitude of the fiuctuations (quantified in terms of as), 
the form of the power spectrum P{k), and the PDF Paly) 
(renormalized to have mean zero and standard deviation 
one) as a function of scale. For our non-Gaussian models we 
do not have simple analytic expressions for these quantities, 
so instead we must measure them directly in the simulations. 
To do this, we output the particle positions and velocities at 
an early redshift zi — {zi — l)/2, when all scales of interest 
are still in the linear regime. We compute the PDF of the 
field on a range of scales R by throwing down 30,000 spheres 
at random and counting the distribution of particle numbers 
in each sphere. We compute the power spectrum by resam- 
pling the particle distribution as a density field on a 128^ 
lattice, and performing the appropriate Fourier transform. 
With Gaussian simulations, it is not normally necessary 



to carry out these steps, since, in constructing a realization 
of a Gaussian field it is straightforward to ensure that it 
has the correct normalization, power spectrum and PDF. 
In order to test the fairness of the techniques we are using 
to measure these quantities in the non-Gaussian case, we 
carry out an identical analysis on a simulation (which we 
denote GA-1.8) with Gaussian initial conditions and a power 
spectrum P{k) oc in a background cosmology with 

f2m = 0.2 and Ha = 0.8. 

The power spectrum and PDF Pniy) with R = 8h~^ 
Mpc (the Lagrangian radius of a typical rich cluster) of 
the linear initial conditions for each model are shown in 
Figs. |l|, 1^, ^, ^, and |^. The solid lines show the quantities 
computed from the simulations, and the dotted lines show 
simple analytic fits that we have used as input to the PS pre- 
dictions. We fit the power spectrum with a two parameter 
model: 



(24) 



where we compute the parameters Ci and C2 for each model 
by a simple least squares fitting to all points k < /cNy, and 
/cNy = nL/up is the Nyquist frequency of the box. The en- 
hanced power in modes with k > fcNy is an artifact of the 
finite mass resolution of the simulations; it does not affect 
the subsequent evolution. For the cosmic string case, use of 
two parameters significantly improves the fit to the data, 
but for all other models we can obtain a good fit by setting 
C2 = 0. In the ICDM-2.4 model, the best fit power spectrum 
has a slope nip — —1.6, which should be compared with 
value nip = 3 -f 2m^ = —1.8 which we would expect from 
a theoretical calculation (see eq. The slightly shallower 
slope arises from missing power on large scales, where the 
assumption of a pure power law in the scalar field (f) breaks 
down due to finite size effects. In the GA-1.8 model, the best 
fit slope is P{k) oc fc~^'*^, very close to the theoretical value 
of -1.8 which we have fed in. 

We fit the PDF PR{y) with a log-normal distribution, 
translated and normalized to have mean zero and standard 
deviation one, that is: 



P-LN(y) = ^ e ' ' 



where 



x{y) = 



\/27rA2 



\n{Cy\A\/A + B) 
\A\ 



with 
B = 

C = y/B^-B2. 



0-4/2 



(25) 

(26) 

(27) 
(28) 



The log-normal distribution has one free parameter A, with 
A > giving rise to an extended tail of positive fluctuations, 
A < 1 giving rise to a suppressed tail of positive fluctuations, 
and Gaussianity in the limit A ^ 0. We fit the PDF Pniy) 
in each model by finding the value of A for which the number 
of peaks of height 3-a or greater is equal to that for Pniy), 
that is we solve for A satisfying 

PR{y)dy= / PMy)<iy. (29) 

/3 J3 

The log-normal PDFs fitted by this procedure also give a 
good fit for all values of y, as demonstrated by our figures. 
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P{k) 






PR{y) 




Model 








o 


Ci 


C2 




A2 




GA-1.8 


0.2 


0.8 


19 


0.98 


-1.82 














ICDM-2.4 


0.2 


0.8 


39 


0.99 


-1.60 





0.58 








ICDM-2.0 


1.0 


0.0 


29 


1.0 


-0.98 





0.35 








Strings-HDM 


1.0 


0.0 


19 


0.58 


-6.95 - 


3.67 


0.17 








Strings-CDM 


1.0 


0.0 


29 


0.86 


-1.9 


-0.6 


0.19 


-0.012 





Voids 


1.0 


0.0 


15 


0.89 


-2.07 





0.154 


-0.049 


0.001 



Table 1. Summary of each of the models, including background cosmology, normal- 
ization, power spectrum and PDF parameters. 





Figure 1. Measured (solid lines) and fitted (dotted lines) power 
spectra (top panel) and PDF (bottom panel) in the Gaussian 
GA-1.8 model. 



Figure 2. Measured (solid lines) and fitted (dotted lines) power 
spectra (top panel) and PDF (bottom panel) in the ICDM-2.4 
model. 



We find that we can fit the scale dependence of the param- 
eter A using a quadratic function in K, truncated to ensure 
that it does not change sign, that is 



A{R) = 



Ax -f A2R + AzBf- 




■ (Ai + A2-B+A3fl^)B=s (30) 

. otherwise. 



For all but the Voids model, PR{y) is very nearly scale inde- 
pendent over the range of scales relevant to cluster formation 
in our simulation, and it is sufficient to set the parameters 
A2 and A3 to zero. In the Gaussian case, we find that the 
best fit value of A is very slightly greater than zero, but that 
the predictions for N^m are virtually identical whether we 



use this best fit value or the choice A — 0. Values for the 
parameters ag, Ci, C2, Ai, A2, and A3 used to fit the nor- 
malization, power spectrum, and PDF in each of our models 
are given in Table Rl. 

In Figs. 0, |, |; 0, 0, and |l| we compare the PS pre- 
diction A*'>A/ for the mass function (curves) with the mass 
function N^J^.J measured in the simulations (data points), 
for a range of output times. For the simulations, we restrict 
our analysis to that portion of the mass function Nyj^f for 
which M > 16Mp, where Mp is the mass of a single parti- 
cle. In the Gaussian case, the PS prediction typically gives 
a good fit to the mass function observed in simulations for 
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Figure 3. Measured (solid lines) and fitted (dotted lines) power 
spectra (top panel) and PDF (bottom panel) in the ICDM-2.0 
model. 



AI > Mnl, where A^nl is the typical nonlinear mass satisfy- 
ing 5c/o"fl(Af„L) ~ ^- show regions for which M < Mnl 
in our figures by using dotted lines to continue the A^>m 
curves. 

A rough estimate of the type of errors we might ex- 
pect in the PS prediction for the non-Gaussian mass func- 
tion comes from consideration of the correction factor / (see 
equation For Gaussian fluctuations, the choice / = 2 
gives a good fit to the number abundance of clusters ob- 
served in N-body simulations, and is justified analytically 
by the excursion set derivation of the PS formula (Bond 
et al. 1991). For the non-Gaussian case, we have again fixed 
/ by ensuring that the mass function accounts for the entire 
mass of the universe, thus maintaining the spirit of the PS 
derivation. However, we cannot be sure a priori that this 
this correction will still account correctly for the number 
abundance of the rarest peaks. For the ICDM-2.4 model, 
the correction factor corresponding to the PDF used in our 
calculations is / = 2.6, suggesting that errors of order 30% 
in the predicted mass function would not be surprising. 

We quantify the goodness of fit of the predictions at 
each redshift by computing the mean and the rms error, 
where the mean error E is given by 

E^-YiN^M^-N^Mj, (31) 

J ^ 

I— 1 



Figure 4. Measured (solid lines) and fitted (dotted lines) power 
spectra (top panel) and PDF (bottom panel) in the Strings- HDM 
model. 

and the rms error V & is given by 

Here Mi is the mass corresponding to the i*'^ data point 
and j is the number of data points. Together, these two 
measurements allow us to quantify both the systematic error 
in our predictions, as well as the scatter about them. Errors 
for each redshift are quoted in Table ^ In computing these 
errors we do not attempt to fit for points for which M < 
i\fNL, since the PS formalism is known to break down for 
these masses in the Gaussian case. We see that typical errors 
are less than of order 0.1 in log^, (Af>M/(ft~^Mpc)^). The 
worst systematic error is for the Voids model at z = 1.7, 
where E — 0.11. The worst rm s e rror is for the Strings- 
HDM model at z = 0, where \^ = 0.18. Even in these 
worst cases, the errors are of order or smaller than typical 
measurement errors in A'^> m , and many orders of magnitude 
smaller than the amount by which predictions for different 
models can differ. 



6 CONCLUSIONS 

We have carried out N-body simulations of a number of mod- 
els with non-Gaussian initial conditions, and verified that a 
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Figure 5. Measured (solid lines) and fitted (dotted lines) power 
spectra (top panel) and PDF (bottom panel) in the Strings-CDM 
model. 
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Figure 6. Measured (solid lines) and fitted (dotted lines) power 
spectra (top panel) and PDF (bottom panel) in the Voids model. 



Model 


z 


E 




GA-1.8 


0.0 


0.006 


0.044 




0.66 


0.080 


0.10 




1.5 


-0.052 


0.097 


ICDM-2.4 


0.0 


0.10 


0.12 




3.0 


0.049 


0.064 




5.7 


-0.051 


0.16 


ICDM-2.0 


0.0 


0.14 


0.90 




1.0 


0.066 


0.11 




2.0 


-0.056 


0.12 


Strings-HDM 


0.0 


0.033 


0.18 




0.33 


-0.052 


0.11 




1.0 


-0.040 


0.16 


Strings-CDM 


0.0 


0.10 


0.12 




1.0 


0.081 


0.10 




2.0 


0.005 


0.065 


Voids 


0.0 


0.047 


0.13 




1.0 


0.057 


0.084 




1.7 


0.11 


0.13 



Table 2. Mean and rrns errors in the predicted mass functions 
at each redshift. 



modified version of tlie Press-Scliechter formalism is able 
to give a good fit to the observed evolution of the cluster 
number abundance over a wide range of redshifts. While 
we have only tested the fit for a finite number of models, 
the models span a range of possible types of non-Gaussian 
behaviour. The ICDM and string models have a skew posi- 
tive PDF, and typically give rise to more clusters of a given 
mass than would be expected in the corresponding Gaussian 
case. The Voids model has a skew negative PDF, and typi- 
cally gives rise to fewer clusters than would be expected in 
the corresponding Gaussian case. The ICDM model is gen- 
erated purely by a local transformation on a Gaussian field, 
so might be thought most likely to agree with the PS pre- 
diction. Tiic Voids model, on the other hand, is generated 
by a more complex non-local transformation of a Gaussian 
random field, and the strings models are generated by an 
entirely different process altogether. 

Quantitatively, the fits to the predicted mass function 
are slightly worse in the non-Gaussian models than in the 
Gaussian one. Typical fits in the non-Gaussian models have 
rms errors in logjQ(Ar>M/(/i^^Mpc)^) of V^E^ ~ 0.1, corre- 
sponding to an error of roughly 25% in N^m- The worst fit 
is in the z = output of a strings model, where the rms er- 
ror is V^B^ = 0.18, or roughly 50% in N>m- These errors are 
of order or smaller than typical observational uncertainties 
in the determination of the mass function, and considerably 



© 0000 RAS, MNRAS 000, 000-000 



Cluster evolution in non- Gaussian models 9 




Figure 7. PS prediction for tiie mass function (curves) and that 
measured in simulations (data points) at various redshifts in the 
Gaussian GA-1.8 model. 



Figure 8. PS prediction for the mass function (curves) and that 
measured in simulations (data points) at various redshifts in the 
ICDM-2.4 model. 



smaller than the amount by which predictions for different 
models may differ. 

Given uncertainties at this level, our results make it pos- 
sible to use observations of cluster evolution to place strong 
constraints on non-Gaussianity in the primordial universe. 
The constraints from current cluster data are outlined in 
detail in Robinson et aJ. (in preparation). Given an inde- 
pendent measure of the matter density of the universe flm 
(which we can hope to gain in the near future from the 
combination of supornovac and CMB observations), clus- 
ter evolution directly constrains the probability distribution 
function (PDF) of primordial fluctuations in the universe. 
If we assume fim = 1, Robinson et al. show that current 
data detects non-Gaussianity, and constrains the primordial 
PDF to have at least 2.5 times as may 3-a peaks as a Gaus- 
sian distribution. If we assume Om = 0.3, current data is 
consistent with Gaussian fluctuations, with the PDF con- 
strained to have no more than 9 times as may 3-cr peaks as 
a Gaussian distribution. Taken together with the substantial 
increase of data on high rcdshift clusters wc can expect in the 
near future, the techniques described here can provide pow- 
erful, model independent constraints on non-Gaussianity in 
the primordial universe. 
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Figure 9. PS prediction for the mass function (curves) and tliat 
measured in simulations (data points) at various redshifts in the 
ICDM-2.0 model. 



Figure 10. PS prediction for mass function (curves) and that 
measured in simulations (data points) at various redshifts in the 
Strings-HDM model. 
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Figure 11. PS prediction for mass function (curves) and that 
measured in simulations (data points) at various redshifts in the 
Strings-CDM model. 



Figure 12. PS prediction for the mass function (curves) and that 
measured in simulations (data points) at various redshifts in the 
Voids model. 
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